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Abstract 

We consider the relativistic Euler equations governing spherically symmetric, perfect fluid 
flows on the outer domain of communication of Schwarzschild spacetime, and we introduce 
a version of the finite volume method which is formulated geometrically (without choosing 
coordinates a priori) and is well-balanced, in the sense that it preserves steady solutions to 
the Euler equations on the curved geometry under consideration. In order to formulate our 
method, we first derive a closed formula describing all steady and spherically symmetric 
solutions to the Euler equations posed on Schwarzschild spacetime. Second, we describe a 
finite volume method which is formulated geometrically from the family of steady solutions 
to the Euler system. Our scheme is second-order accurate and, as required, preserves the 
family of steady solutions at the discrete level. Numerical experiments are presented which 
demonstrate the efficiency and robustness of the proposed method even for solutions containing 
shock waves and nonlinear interacting wave patterns. As an application, we investigate the 
late-time asymptotics of perturbed steady solutions and demonstrate its convergence for late 
time toward another steady solution, taking the overall effect of the perturbation into account. 
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1 Introduction 

The finite volume method is a versatile technique for scientific computing, which has found many 
applications in physical and engineering sciences. In particular, it allows one to approximate weak 
solutions (containing shock waves) to nonlinear hyperbolic systems of balance laws such as, for 
instance, the Euler equations of compressible fluid dynamics. In the present paper, we propose 
a geometric version of the finite volume method for general balance laws of hyperbolic partial 
differential equations, and we apply this method to the Euler equations for spherically symmetric, 
relativistic fluid flows posed on a curved spacetime and, for definiteness, the outer domain of 
commimication of Schwarzschild spacetime. The proposed method is second-order accurate (in 
smooth regions of the flow), and well-balanced in the sense that steady solutions to the relativistic 
fluid equations are preserved at the discrete level. 

The balance laws of interest in the present work have the following general form. Given a 
spacetime (M, g) with Lorentzian metric g and covariant derivative operator V, we consider the 
class of balance laws 

ya(T"p{(p)) = 0, (1.1) 

where T"p{(p) represents the energy-momentum tensor of an (unknown) tensor field cp defined 
on M (the indices a,jS ranging between and 3). We use a standard notation for the metric 
g = gafidx"dxf^ in coordinates (x"), and repeated indices are implicity summed up. We lower (or 
raise) indices with the metric gap (or its inverse g"^^) so that, for instance, u„ = gaput^ for a vector 
field u". 

In particular, we are interested in the relativistic Euler equations for perfect compressible 
fluids, corresponding tocp = (p, u) in | |1.1) with 

T"p{p,u) = {pc'+p)u"up+pg"p. (1.2) 

Here, the scalar field p > denotes the mass-energy density of the fluid and the vector field 
u" its velocity, normalized so that u"u„ - -1, while c > represents the light speed. Moreover, 
the pressure function in {l.2\ is given by an equation of state p - p{p) which must satisfy the 
(hyperbolicity) condition p'(p) e (0, c) (for all p > 0), so that the equations | |1.1) can be written in 
local coordinates as a system of nonlinear balance laws, which is strictly hyperbolic for p > 0. In 
general, initially smooth solutions to |[LTJ-|[1]2} become discontinuous in finite time and shock 
waves form and then propagate within the spacetime. 

A broad literature is available on the design of robust and accurate, shock-capturing schemes 
for general hyperbolic systems posed on a flat geometry like the Minkowski spacetime. In the 
present work, we intend to also take a curved background geometry into account, by following 
recent work by the first author and his collaborators; cf. |[2l [3] [T2|. To this end, we introduce 
a finite volume scheme which is based on the geometric formulation | |1.1[ |, rather than on the 
corresponding partial differential equations in a specific local coordinate chart. Moreover, in 
order to achieve the well-balanced property, we extend the approach in Russo et al. EOl |2T1 l22l 
and LeFloch et al. [11 j, and we introduce a discretization which accurately takes into account the 
family of steady solutions to the balance laws and, therefore, the geometric effects induced by 
the Lorentzian geometry (M, g). To implement this strategy, it is necessary to first investigate the 
class of steady solutions to the Euler system on the curved background under consideration. 

Numerical relativity has undergone a tremendous development in recent years, and the reader 
is referred to Marti and Miiller [,151 for a review of numerical methods developed first for special 
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relativity, and to Font ||7| and Alcubierre ||T| for a review in the context of general relativity. 
See also Papadopoulos and Font flS^ 191 and Novak and Ibahez flT'l. Various astrophysical 
applications have been successfully dealt with in recent years, including the evolution of neutron 
stars and the merging of black holes. For background on general relativity, we refer the reader to 
p3] and, for theoretical and numerical tools concerning nonlinear hyperbolic equations and their 
discretization, we refer to [4, 9 J and [14|, respectively. 

An outline of this paper is as follows. In Section|2j we introduce a version of the finite volume 
method, which we state fully geometrically on a curved spacetime and for the general balance laws 
In Section |3j we consider the relativistic Euler equations when the background is chosen to 
be the outer domain of communication of Schwarzschild spacetime, and we determine all steady 
solutions to | |1.1| -( [T!2) in this context. Next, in Section]!] we introduce a finite volume method for 
the Euler system, which is well-balanced and second-order accurate. Finally, in Section|5| various 
numerical tests are presented, which demonstrate the efficiency of the proposed scheme; as an 
application, we study the late-time asymptotics of steady solutions under perturbation. Section|6] 
contains concluding remarks and refers the reader to the follow-up paper LIOJ concerned with 
self -gravitating matter in spherical symmetry. 



2 The geometric finite volume method on a curved spacetime 
2.1 Spacetime foliations by spacelike hypersurfaces 

We begin by presenting a framework based on a general four-dimensional spacetime (M, g), 
that is, a manifold (possibly with boundary) endowed with a Lorentzian metric g with signature 
(-, +, +, +). We denote by V the spacetime Levi-Civita connection associated with this metric. As 
is customary, we assume that M admits a foliation {^^fj^^jp ^ by oriented spacelike hypersurfaces 

such that the parameter t : M. —> R+ provides us with a global time function, satisfying dt i= 
on M. This allows us to distinguish between future-oriented {t increasing) and past-oriented 
{t decreasing) timelike directions on M. By definition, we thus have M = Uf>o ^t, and M is 
topologically diffeomorphic to 1R+ X IKq while its boundary is the imion of the initial slice IKq and 
a boimdary R+ x d^o determined from d^o, i.e. 

M U = iKf, 5{f = Jft U d:Kt. 

t>o 

(Observe that the spatial slices may have a non-trivial boundary.) 

Given local coordinates {x") = (f,x^), in which Greek indices describe 0, ...,3 while Latin 
indices describe 1, 2, 3, we express the spacetime metric in the form g = gap dx"dx^, and we denote 
by {g"l^) its inverse. The assumed (3 + l)-decomposition of the spacetime is standard in general 
relativity (cf., for instance, the textbooks 11, ,23,1 ) and is determined by the time function t. We 
also choose coordinates (xi) on the initial slice JCq and propagate them in the spacetime along the 
vector field Vf. This leads us to the metric decomposition 

ds^ = -N^df- + g, (2.1) 

where N = {-g{dt, dt))~^^-^ > is referred to as the lapse function and 'g = 'gt = gij dx'dxi represents 
the induced Riemannian metric on the slices. Denoting by dVg and dVg^ the volume forms 
associated with the Lorentzian and Riemannian metrics, respectively, we can write 

dVg=NdtdVg, B = [detgf'^, NB = {detg)^'^. (2.2) 

Let n = NVthe the future directed, timelike unit normal to the slices and K be the second 
fundamental form, defined by K{X, Y) = -g( Vxw, Y) for all vectors X, Y tangent to the hypersurf ace 
Mf, so that 
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in local coordinates {t,xi) adapted to the foliation. We introduce also the Levi-Civita connection 
V of the slices {Dit, gt), given (for any tangent vector fields X, Y) by 

VyX = VyX + K(X,Y)M. 

Finally, the Christoffel symbols F^^ = j g^'^l^agpe + d^gae - dggap^ of the spacetime metric read 

r° - — F^' - -Q'''dM^ 

1 1 dB (2-^^ 

2.2 Spacetime formulation of balance laws 

Consider now the general system of balance laws and, by first assuming enough regularity 
on the solutions, let us rewrite it in local coordinates adapted to the (3 + l)-decomposition 
ijZTJ determined by the time function f. From the definition of the covariant derivative, is 
equivalent to 

d^T^?' + d-T^?' + FgoT"/' + FjoT'"' + Yip?' + FJ.T'^'^ 

or, in view of the expressions of the Christoffel symbols \2.3\ , 

doT"° + djV" = -da{ln{BN^)) T*' - do{lnN) T"" + T'i, 
doT°' + djV' = -da{ln{BN^)) V" + ^g'''dk{N^) T"" 

+ 2N - F' V'' + 2 ^ T"' + 2 ^ TL 

' I'' N B 

After multiplication by the weight BN^, we obtain the following formulation of the balance laws 

do{BN^T°°) + dj{BN^T°') = S°, 

doiBN^T'^) + djiBN^ri) = S\ ^^'^^ 

with right-hand sides 

S'=BN^d,NT''>-^d,g,jT'^K 

S' = ^g"dkN^ TOO _ gj^3 ^,kg^^,^ jjo _ gj^3 ^i^jjk 

- T°'do{BN^) + 2T°'BN^doN + IN^T'idjB. 

In particular, plugging in 1 2.4 1 the expression of the matter tensor | |1.2[ | yields the formulation of 
the Euler equations on a curved spacetime: 

do(BN^({pc^ + p) m° - ^p)) + <9;(B]V3 (pc2 + p) u° ui) = S", 
dofSN^ipc^ + p) M° m') + dj(BN^({pc^ + p) u' ui + p g'')) = S', 

which consist of four equations for the five unknowns (p, u"), satisfying the constraint u"Ua - -1. 

Introducing local coordinates {x") = {t, xi) and recalling the decomposition | |2.2) of the volume 
dVg, the formulation dZSl can be recovered in the sense of distributions. 
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2.3 The geometric spacetime finite volume method 

We are now in a position to introduce the geometric formulation of the finite volume method, 
which we design directly from the covariant form of the balance laws, rather than introducing 
first coordinates and then a discretization. 

First of all, we introduce a triangulation of the spacetime, say 

m = [Jk, 

made of finitely many open sets K, which are assumed to satisfy the following conditions: 

• The boundary dK - Ucc^k^ is piecewise smooth and contains two spacelike faces (i.e., 
having an induced metric of Riemannian type) denoted by and e^, and timelike (or 
"vertical") faces (i.e. having an induced metric of Lorentzian type), the latter faces being 
denoted by 

e° ed°K = dK\[e^,e-]. 

• The intersection K n K of two distinct elements is a common face of K and K', or else is a 
smooth sub manifold with dimension at most 2. 

We then adopt the following notation: 

• Along the timelike faces e^, we introduce the outgoing unit normal vector field denoted by 

flK- 

• |e^|, |e~|, |e°| denote the Lebesgue measure of the sets X,e^,e~,e°, respectively, which is 
defined from the Lorentzian metric or the induced metric on these hjrpersurfaces. 

Furthermore, when the spacetime is endowed with a foliation by spacelike hypersurfaces, 
say M U (9M = \Jt>o ^t, associated with a time function f : M IJ dM [0, +oo) we say that the 
triangulation 

MudM=[jK, 

is compatible with the time function t if it is determined from a sequence of discrete times 
and a triangulation T' of the initial three-dimensional slice "Kq, say 

in such a way that the boundaries of the elements K are transported to the whole spacetime along 
the vector field dt associated with the time function so that all the vertical faces are parallel to this 
vector field. This property makes it clearer to advance the numerical solution forward in time, 
and is assumed from now on. 

The finite volume method is then based on the following general weak form of the system of 
balance laws. Recall that solutions to (nonlinear hyperbolic) balance laws are generally discon- 
tinuous, and these equations must be understood in the sense of distributionals. Hence, we seek 
here for weak solutions for which is understood in the averaged sense 

r 7i|Jr''rfy, = o, (2.6) 
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in which X" is a test-field (i.e. is smooth and compactly supported in M) and tt'^^ denotes its 
deformation tensor defined by 

nfi! = l{y.Xp+ypX,,). (2.7) 

The finite volume method is based on the above integral formula, except that we must now 
integrate over an arbitrary spacetime element K e 7. Given any smooth vector field X, which no 
longer needs to be compactly supported, we write 

£Tix"nK,pdV,. 

= [ TU"nK,pdV,.- y [ Tix"nK.pdV,.^ + f nfh^S'dVg, 

in which obvious notation has been used for the induced volume form along each boundary 
component of K. The approximation scheme is now defined from this general identity, by 
choosing X to be either the vector dt associated with the time function, or vector fields tangent to 
the spacelike hypersurfaces. 

Under our assumptions, the normal n^.p along the spacelike sides has components (-N, 0, 0, 0), 
hence the above equation becomes 



X 



r TO'^X^NdV^- Y r Tix"nK,pdV,o^+ f nfh"^dVg. 



Finally, in specifically chosen coordinates, we can choose covector fields with constant compo- 
nents, say X*"' = (1, 0, . . .), X'^' = (0, 1, 0, . . .), etc., and we can introduce the source-terms 

which allows us to express the following averaged balance laws (for a = 0, 1, . . .) 

1 T°"NdV,.= f T°"NdV,.- Y f T"l^ nK,pdV,o^ + f S" dVg. (2.8) 

The geometric version of the finite volume method is based on the integral identity | |2.8| . The 
solution is represented, at every discrete time t„ and within each spacelike hypersurface (say, ?{(„) 

by constant states T^- = [T^-^ y The constant states Tp+ on the "next" h5^ersurface (that is, iKf^^J 

are then determined as follows. Along each vertical face of an element K, we a priori fix a 

family of numerical flux functions, say, F"^ [T^-^ , T^^ j (with K' defined by K n K' = and K + K'), 

for the approximation of the vertical contribution, so that 

tn,,p^P^,^(f:i,f,l). 

The numerical fluxes are assumed to be locally Lipschitz continuous and satisfy standard proper- 
ties of consistency and conservation ||3l. The finite volume scheme generates the constants values 
— o« 

Tg+ and reads (for a - 0,1,...) 
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in which Nc± are consistent approximations of the lapse function on the corresponding hjrpersur- 

—a 

faces and S are consistent approximations of the source-terms. For instance, under the symmetry 
assumption of main interest in this paper, we can work in the quotient manifold which has space- 

—0 „ —00 „ —01 

time dimension (1 + 1), and we need to introduce consistent approximations of S - ^qqT +^^^7 

-1 1 —01 -, —11 

andS =rJjT +T\J . 

Furthermore, a restriction on the time step is required for stability purposes. An extension 
of this scheme will be introduced in Section 4 below in order to make it to preserve steady state 
solutions at the discrete level. 



3 Steady fluid flows on Schwarzschild spacetime 
3.1 Steady solutions on a curved spacetime 

Our first task is to investigate the properties of steady state solutions to the Euler equations \2.5\ . 
Without imposing symmetry assumptions, it is clear that no analytical closed formula could be 
derived for the solutions to 

dj(BN^ {pc^ + p) u° u') = S°, 

(3.1) 

dj{BN^[{pc^ + p) u' u' +pg'')) = S', 

which is a mixed type (hyperbolic, elliptic) system in the variables (x^). By imposing spherical 
symmetry (for instance), this system reduces to a system of two ordinary differential equations, 
which, in itself is already quite chalenging. 

di(BN^ ipc^+p) u°u^) = S°, 

(3.2) 

d,(BN\{pc^+p){u'f+pg'')) = S\ 

Furthermore, we point out that the solutions studied in the present section were first introduced 
in the physics literature in a different gauge [8J and represent the steady state accretion of matter 
in a Schwarzschild black hole geometry. 



3.2 Euler equations on Schwarzschild spacetime 

We are now interested in the outer domain of communication of Schwarzschild spacetime, which 
describes the exterior of a spherically symmetric black hole. The Schwarzschild metric is a 
particular solution to the Einstein equations and, in the so-called Bondi coordinates (f, r, 6, (p), 
reads 

= -(l - de + (l - ^)"' dr^ + r\de^ + sin2 (3.3) 

which is meaningful for r > 2m. The coefficient m represents the mass of a black hole located at 
r = 0. This spacetime is spherically symmetric, that is, is invariant under the group of rotations 
acting on the spacelike 2-spheres of constant t and r. It is static, since the vector field = 
is a timelike Killing vector and asymptotically converges to the (asymptotically flat) Minkowski 
spacetime, when r +oo. The expression | |3.3^ represents the outer domain of communication, 
only, and a different coordinate choice would be required to go through the horizon located the 
"boundary of the coordinates" r = 2m (around which the spacetime itself is actually regular). 
From now on, for simplicity in the presentation we assume that the equation of state is linear, 

i.e. 

Pip) = a'p, (3.4) 

where a is a constant. This is not an essential assumption, however. We restrict attention to 
solutions depending on the radial variable r, only, and such that the non-radial component of the 
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velocity vanishes. In orther words, we have = (u^{t,r),u^{t,r), 0,0^ and, as a consequence, 
the energy momentum tensor satisfies T^^ = T"^ = T^^ = T^^ = T^^ = 0. 

The velocity m is a unit vector, thus -1 = -(l - ^) (m°)^ "*"(■"■" ^) ("^^^ ^ terms of the 
rescaled velocity component V = ^zs: ^, we find 



(m0)2 ^ C_ ^^1^2 ^ ^2 
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We can thus express the Euler system in the form 

do{BN^T°°) + dr{BN^T^°) = S°, 
do{BN^T°^) + driBN^T^^) = S\ 

in which 

S° = 0, BN^ = sine rir- 2m), 

gi ^ joo^^n^^^i _ BN^V' + 2N'T'%B. 

We thus arrive at 

dt (r(r - 2m)P°) + dr {r{r - 2m)T°'^) = 0, 

2 

dt (r(r - 2m)T°i) + dr (r{r - 2m)r") - 3mT" + ^ (r - 2m)^ 1°° 
- r (r - 2m)^ T^^ - r sin^(0) (r - 2m)^ T^^ = 0, 

whose coefficients are given by 

?flo _ /i _ ?^Woo _ g'p + P(P)^'/g^ 2 ^01 _ 7^01 _ cV + P(P) 

_ 2m J c2 - y2 ' 



and 



22 _ pip) t33 _ P(P) 



r2 ' r2 sin^e' 

In conclusion, the Eu]er system on a Schwarzschild background takes the form 

dt (!l!:^?oiJ ^ _ 2^)2^11) _ 3^ 'l^^u (3.6) 

+ m 1°° (r - 2mf r- = 0. 

r r - 

For later use, we record here some additional formulas: 



l + ff2y2' (3.7) 
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and 

jOQ jU _ _ ^2 ^2 p2 jQO _|_ jU _ (g + C ) + (V + C ) 

4_v2"2^' ' (3.8) 

TOO _ jU ^ (,2 _ ^2) (^11)2 _ (^01)2 ^ ^-^-^ p2^2_ 

3.3 A closed formula for steady fluid solutions 

We now derive a closed formula for smooth steady solutions p = p{r) and V = V{r) to the Euler 
equations posed on a Schwarzschild spacetime. Various plots of solutions are provided in Figures 
3.1 to 3.4. We work in an interval r e {2m, R] for some fixed R > 2m, and we impose boundary 
data at r = K 

p{R)>0, V{R)e{-l,l). 
The Euler equations reduce to the following system of ordinary differential equations: 



d (r{r- 2m) -^^ 



dr 



T"' = 0, (3.9) 



d , _ , _ _ jOO _ Til 

^ (r{r - 2m)T" - otT" + mT°' - 2a\r - 2m) , , = 0. (3.10) 
The "first" equation ||3.9^, after integration over the interval [r, R\, yields 



c^p + vipir)) 

A{r) = r{r-2m) ; \)y'' V{r)=A{R), re{2m,R), (3.11) 

c2 _ y(j-)2 

where A{R) is determined by boundary data prescribed at r = R. This implies that if V{R) % 0, 
then V{r) § for all r e {2m, R). By solving | |3.1H in terms of V{r), we deduce 

-KP+ J{k p)2 + 4c2A2 
V = ^ , (3.12) 

where k = K{r) - r{r - 2m){c^ + o^) is a function determined by the mass and the sound speed. 
By using our expressions p.8^ , the equation | |3.9^ (after taking the square) is equivalent to 



/ r^r-2m)^ /~oo _ o^c^ -qq _ 2\\ ^ 
dr\ c2 I (c2-a2)2^ Wj 



Hence, by introducing the following suitably weighted quantities 

T»o = r(r-2m)T™ T" = r(r - 2m)T", e=^^, a=—, 
the Euler system becomes 

d_ ^ JLpo T" - e2(T00 - T")2)) = 0, (3.13) 

i.Tii + 1 f .JH^ _ . (TO" - T") = 0. (3.14) 

We can integrate the equation 1 3.13) , as we did earlier, and with the new notation we now 
have 

T00Tii-e2(T00_^ii-,2^^2^2^ 

where we recall that A = A{R) = A{r) is a constant. By setting Y = T^^, this equation takes the 
form ^ ^ 

-e2(T™-Y)2 + T00Y-c2A2=0, 
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which leads us to 

T"" = Y + l) + ^ Vy2(l + 4^2) - ieVA^. (3.15) 

In addition, since the dominant energy condition (satisfied by the fluids under consideration in 
134)) imposes > T", we must have f°° -Y> 0, that is, 

which is a constraint on the values taken by Y. Using the definition of Y and the expression 1 3.15^ 
in the "second" Euler equation | 3.14^ , we now find 



dY 1 I m 
dr r\r- 2m 



(7^~^'')M^^ Vy^(l+4e2)-4e2c2A2) = 0. (3.16) 



This differential equation in Y can be solved and we obtain 

dY Urn 



Y + VY2(1 + 4e2) - 4e2c2A2 r\r-2m 
and, after integration. 



^e)^^dr 



J, Y + VY2(1 + 4e2) - 4e2c2A2 ~ 4^2 J, I r {r-2m)j 

-"((7)^(^)1 



F{R)- 




100 1000 10000 



T-{11) 

Figure 3.1. Steady solutions for three values of the radius r. 



With the change of variable -Y^A^ = Y^ ^J^^fgi , the integrand simplifies: 

d(coshX(y)) _ 1 r« e2X(y) _ ^ 
" J, cosh X(y) + b sinh X(y) " (1 + 6) J, e2X(y) + C 
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Figure 3.2. Two steady solutions Y = T^^ for the radius r = .8 
with p{R) = 460 and ViR) = .001. 
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Figure 3.3. Two steady solutions T^^ for a fixed radius r = .8 
with p{R) = 460 and V{R) = .001. 
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where we have used the change of variable y = cosh X{y) and set 6 = (1 + 4e^)^^^ and C = jif • We 
thus obtain 

with e^xCy) _ ^^^^ integration, 

= ^(x(y(r))-X(y(R))) 

1 / / X e2X(y(R)) + r 

+ 2(x(y(r)) - X{y{R))) + In 



e2X(y(f)) + Q 

/I 2 \/ \ 1 / e^^W + 

= ( - ^ T35^)K)) - W))) + ^ m ^.j^,;^ ^. 

^ 1+0 ' 

Therefore, we have derived an algebraic relation for the function X = X(j/(r)), that is, 
= 0, 

where we recall that 

^_ OC , , 7,1/9 ^ 1+6 ^, , /1\^ 



y(R)+VyW^/ ly(r)+V^(^+i|§/ If(R)/ (3.17) 



T 1/T 1+6 /1\ \* 

^, 6 = (l+4e2)i/2, C=^, ^^''^ = (r) V-^mj 



In summary, given p(K) > and e (-1, 1), there exist two steady solutions 

y± = y±{r), r e (2m, K] 



12 



to the reduced Euler equation | |3.17[ |. We have derived an algebraic relation for these solutions, 
specifically H{y{r)) - 0, so that y{r) can be computed by solving this equation for each r, for 
instance by a fixed point technique. From the function y = y{r), it is then straightforward to 
recover the physical variables. Namely, by using | 3.15[ , 



= y + l) + ^ Vy2(l + 4^2) - 4e2c2A2, 



and T^" - T^^ = {c^ - o^) p, we obtain the expression of the density p = p{r) and, next, for the 
velocity V = V{r) from \3.12\ , that is. 



V = 



-Kp+ yjK^ p2 + 4c2A(K)2 

2A{R) ■ 



4 Well-balanced approximations on Schwarzschild spacetime 
4.1 Finite volume method 

We express the Euler equations on Schwarzschild spacetime in the form of a hyperbolic system 
of balance laws, that is, 

dtV + (9,.F(U, r) = S(U, r), r > 2m, (4.1) 



with, in view of l|3.5b-l 3.6 1, 



U = 



9 cV+p(p)VVc^ \ 



_ 2m) V 



and 



F(U,r) = 



S°(U,r) 



F''(U, r)\ _ ( r(r - 2m) 



FHV,r)j [(r-2m) 
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^2 V>+p(p) 2 



p ((3ot(c2 + 02) - 2a2r) y2 _ c^fmic^ + o^) - 2o^r)) 



~ I r-2m 
r(c^-V^) 

We apply the finite volume technique, presented in Section [Z3] for general systems of balance 
laws, by working here with the (1 + l)-dimensional quotient (by the group of spatial symmetries) 
of the Schwarzschild spacetime. For this quotient metric g, we thus have ds-^ = -c2 ^1 - dfi + g 
with induced volume form 



\l/2 



dr = [1 j dr, 



-1/2 



dr. 



In agreement with Section 2.3 and by introducing the approximations 



and 



—n 1 r>*'l^ -n 1 r'> 

r'*'''/ 2m X -1/2 

= AV^„ = (i - -) dr 

-''•,-1/2 



S(i„,r)rfyj 



the finite volume scheme takes the form 



— n+l — n l^t i—n —n \ — n 

=U;--r-(F,,v2-F,-i/2) + AfS,., 



(4.2) 
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in which fj+\/2 are consistent approximations of the exact flux of the system | |4.1) . 

As usual, the Courant-Friedrichs-Lewy (CFL) condition is imposed on the time step in order 
to guarantee stability. Specifically, we set At = t„+i - t„, which we assume to be independent of n 
for simplicity, and impose the inequality 

^ max|A(U)| < 1, (4.3) 

Ar 

where the maximum is taken over the (real) wave speeds A(U) of the Euler system. 

4.2 Taking the Schwarzschild geometry into account 

In steady state solutions to ( |4.1) , which by definition satisfy 

5,F(U,r) = S(U,r), (4.4) 

the source terms exactly balance the flux terms. We will construct the well-balanced version of the 
finite volume scheme introduced in the previous subsection by imposing that the same property 
must hold at the discrete level of approximation for the family of discrete steady states. For 
instance, cell-centered evaluation of the source terms, generally, do not ensure the preservation 
of these discrete steady states. Therefore, we look for an adapted discretization of the source-term 
which directly uses information from the steady state equation and, more specifically, uses the 
characterization ( |3.17| exhibited in Section 3, above. In turn, our scheme will satisfy a discrete 



version of the steady state system 1 4.4 1 



Motivated by the work by Russo et al. EOllZTI , the reconstruction scheme proposed now takes 
the family of steady solutions into account for the numerical evaluation of the intermediate states 
at which the numerical flux is computed. Specifically, recalling the first Euler equations (in the 
form adopted i the present paper) contains no source-term, we define the well-balanced finite 
volume scheme 

— 0,11+1 — 0,n At /— 0,n — 0,fi \ 

— ttI/I i^t /— 1,K — l,n \ . —1," 

^ -^^{FnV2-F,-V2) + ^tSj , 

by substituting our closed expression of the family of steady solutions. Clearly, the first equation 
is in a conservation form and does not need any well-balanced correction, so we concentrate on 
the second equation. 

First of all, we approximate the solution in each cell by a steady solution and, under this 
approximation and by using integration by parts, we can then transform the source term for the 
second Euler equation (cf. the expression | |3.6[ l), as follows (for exact steady solutions): 



\3m T" - m T°° + (r - 2mf — ^ M^tt 



I r 0+1/2 / 

= J- r^dMr-2nrfT^^)(l-'-^)-'''dr 
This identity, valid for exact steady solutions, motivates us to propose the following approxima- 



14 



tion for the source term 

— ll,n — ll,n 



(4.6) 



"^'■;-l/2+ 



Of course, in order to be able to make use of the above definition, it remains, on one hand, to 

— — 11,H 

introduce suitable approximations Tj+112- and T^-i/2+ (consistent with the values taken by the 
"true" solution) at the interfaces between the cells and, on the other hand, to evaluate the integral 
term above. 

Consider first the latter issue, we propose here to use Simpson's rule, and we replace the 
integral term by the following explicit expression 



^ rrj+i/i- — ll,n —n,n —11,11 



(4.7) 



Observe next that, since we have introduced new states at the interfaces between the cells, it 
is natural (and actually necessary in order to achieve the weel-balanced property) to compute the 
numerical flux (of the second Euler equation) in terms of these interfaces value. In other words, 
we write (at j + 1/2, say, and with similar formulas at / - 1 /2) 

fy+l/2 = F {Uj+l/2-' ^+1/2+)' 

-1 _-i/-o -1^-1. (4.8) 

^j+l/2 - ^ y-^j+1/2-' ^i+l/2-' ^j+l/2+' ^j+l/2+}' 

where we are taking in to account the particular dependency of the flux of the Euler system. By 
taking the Schwarzschild geometry into account, we can write (with obvious notation) 

-l,n ^1'" 



^■+1/2- = ry+i/2(ry+i/2 - 2m) T^+1/2-, 
-l,n 

Uy+i/2+ = ry+i/2(ry+i/2 - 2m) T^+i/2+, 



It now remains to compute the states at the interfaces. First of all, recalling that, for exact steady 
solutions, the expression r{r - 2m)T^^ is a constant, we naturally determined the "reconstructed" 

— 01,« — 01,n 

states T 1^1/2+ ^y+i/2+ the interfaces (by interpolation from the states within the cells) by 
setting 

-oi,(! (^rj + Arj){rj + Arj - 2m) 

^^■^1/2. - rj,y2{rj,i/2-2m) ' 

~oi.n ry(ry - 2m) -01." 

'^^■^1/2-- ry^i/2(r;+i/2-2m) ' 

=01." ry(ry - 2m) =oi'« 

^^-1/2^" r,_i/2(r;-i/2-2m) ' 

-oi,H (y. - Arj){rj - Arj - 2m) 

^^-1/2- = rj.y2{rj-v2-2m) ^^"^ ' 

—11,11 — ll,n 

On the other hand, the "reconstruction" of the interface states Ty+j^2+ ^nd Tj^i/2+ is more delicate 
and requires the full algebraic relation ( 3.17[ , which provides us with a complete characterization 
of steady solutions. This characterization is based on the function H = H{y) discovered in 
Section 3 and, therefore, in each computational cell we now impose the two relations (with 
obvious notatopn for the quantities y", y"_i^2+' ^^'^ y/+i/2-) 

Hifj+i/2-) = ^iy") = Hiy"j-i/2+)' 
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in which the value y" is explicitly known from the states in the cell /, while the unknowns i/"_i/2+ 
and y"+i/2- interfaces are obtained by (numerically) solving the above algebraic equations 

(by a fixed point technique). Thanks to the relation | |3.15| , we can now determine all of the 

— 00,« — 00,n 

interface values Tj±ij2+ arid ^y±i/2+- Finally, from the identity | |3.8[ l, we can also compute the mass 
energy density states p"+i/2+ as well as the velocity components Vj^^/2+' and ^^+^^2+- have 
now completed the design of our well-balanced scheme (at the first-order of accuracy). 

In summary, the above construction relying on steady solutions in order to define the interface 
states for the evaluation of the numerical fluxes, it can be checked that steady solutions are exactly 
preserved at the discrete level of approximation. On the other hand, the consistency property of 
the original scheme (ensuring that limits of the scheme do satisfy the Euler equations in the sense 
of distributions) also holds. 



4.3 Second-order accuracy 

To arrive at a second-order scheme, we follow Nessyahu and Tadmor ||T6ll and introduce a 
predictor-corrector scheme, as follows: 

—11+1/2 —n At -' 

=''>-2Ar^r 

— H+l 1 /— (1+1 — H \ 1 ,—' —' , At /—/—n+l/2\ — /— (i+l/2\\ 

U,+i/2 = 2(U; +U,.,0+8(U,-U,+i)-^{f(u,+i )-f(u,. )). 
— / — / 

The states Uy, Fy represent first-order approximations of the space derivatives (of the field and 
flux variables) at the point Vj and can be computed in several ways. A standard choice (which is 
used in this paper) is given by 

— / / — II — n — 11 — 17 \ 

— / /—n —n —n —n \ 

F^. = m(f^.^i-F^.,F^--F^_i), 
where M(U, W) is the min-mod limiter (which we apply composent-wise) 



M(!J,W) = 



I sgn((J) min(|(J|, |W|), sgn(U) = sgn(W), 

1 0, otherwise. 



5 Numerical experiments and applications 

5.1 Comparison between several schemes 

In the following numerical experiments, we investigate the proposed scheme for the computation 
of weak solutions to the Euler equations on Schwarzschild spacetime. We work within the exterior 
domain of communication r e {2m, R) limited by the horizon r = 2m and a sphere with radius 
R > 2m. In all tests, the sound speed is taken to be a = 0.3, the light speed is unit, the upper 
space bound R - 2, and the mass parameter is ?fJ - 0.1. More precisely, our computations take 
place in an interval r e (ro, R), with Tq > 2m, so that we stay away from the horizon, on which the 
Euler equations (in the chosen coordinates) are singular. In every test, we treat the boundaries at 
r = Tq and r = Rhy solving a Riemann problem between the boundary data (determined from the 
given initial data) and the current numerical values at the boundaries, and we use the flux of the 
Riemann solutions, which is the Godunov scheme at the boundary ||5l|6l. 

We begin by illustrating the interest of the well-balanced property and we compare together 
two schemes, the proposed well-balanced one as well as a "naive discretization" (cf . next para- 
graph) of the right-hand sides of the Euler equations. Throughout, we also plot the exact (or 
asymptotic) solution when available. The "naive discretization" of the right-hand sides of the 
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Euler equations is defined by replacing the discretiation of the source Sy in < |4.2[ by an evaluation 

— n — n 

computed from the state \Jj , that is, for the naive version, only, Sj reads 



The numerical results are plotted in Figures 5.1.1 and 5.1.2, where we have chosen an initial 
data made of a perturbation of a steady solution. This solution is determined from the value of 
the velocity V{R) = 0.001 and the density p{R) = 460 at the end point of the interval {ro,R). A sine 
perturbation is added to the steady solution, so that a genuine evolution in time now takes place 
and we work in the interval tq - 0.5 < r < K = 2. We use the space mesh size Ar - 0.025 and a 
CFL number equal to 0.9. 

Figure 5.1.1 (a) represents the steady solution (velocity component) together with its perbur- 
bation which serves as our initial data. Figure 5.1.1 (b) represents the steady solution together 
with the numerical solutions with our two schemes. Observe that the well-balanced scheme pro- 
duces a solution which is closer to the steady solution and oscillates about it, while the standard 
scheme deviates significantly from it. Figure 5.1.1 (c) represents the same solutions, but at a 
much later time: we now observe that the well-balanced numerical solution slightly oscillates 
about the steady solution, while the standard numerical solution is clearly completely wrong. 
Figure 5.1.1 (d) represents the time-asymptotic behavior, and we observe that the well-balanced 
scheme has re-converged to the original steady solution, the perturbations having cancelled out 
asymptotically, while again the standard scheme has generated a completely wrong solution. 

Figure 5.1.2, instead of the physical variables (like the velocity above), shows the nonlinear 
expression r{r - 2m)T'^^, which is known to be constant for steady solutions. Figure 5.1.2 (a) 
represents this function for the steady solution (which is thus a constant) and for the well- 
balanced scheme after N = 500 and N - 1000 iterations, respectively. Figure 5.1.2 (b) is a plot 
of the relative error for the same quantity. Again, we observe that the well-balanced scheme 
produces a quite satisfactory result with 0.5% at the end point r - Rof the interval. The accuracy 
is better near the horizon but grows with r. (A further correction of the scheme may be found 
useful to improve the accuracy for large radius r.) 

5.2 Propagation of a shock/rarefaction pattern 

We study here initial data containing a shock separating two steady solutions. As we evolve this 
initial data, the profile of the solution changes and additional waves arise. The initial jump we 
choose being arbitrary, a full solution to the Riemann problem is generated and both shocks and 
rarefactions may occur. The two steady solutions are defined as follows: the left-hand steady 
solution has the velocity V{R) = 0.001 and the density p{R) = 460, while V{R) = 0.004 and p{R) = 
480 for the right-hand solution. As in the first test, we work in the interval Tq = 0.5 < r < K = 2. 
We use the space mesh size Ar - 0.025 and a CFL number equal to 0.9. 

Figure 5.2 (a) represents the initial discontinuity separating the two steady solutions, while 
the other three plots in Figure 5.2 (b), (c), (d) show the numerical solution given by the standard 
and the well-balanced schemes. We observe that the well-balanced scheme produces a sharper 
solution with a jump of the same magnitude as the initial jump, while the standard scheme has 
generated a "spike" of much larger amplitude. Yet, this unphysical spike is diminishing as time 
evolves and gets back closer to the amplitude of the well-balanced numerical solution. 

5.3 Late-time asymptotic stability of steady solutions 

We have now validated our well-balanced scheme and this motivates us to now apply it in order 
to study the nonlinear stability of a given steady solution. The initial data is defined by adding 
a compactly supported perturbation. The velocity and density at the right-hand point of the 
spatial interval are chosen to be V{R) - 0.001 and p{R) - 460, respectively. We now work in the 
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(a) Initial data and perturbation (velocity) (b) Intermediate time 




(c) Later time (d) Asymptotic behavior 

Figure 5.1.1. Numerical solutions with the standard and well-balanced schemes. 



interval tq = 1.2 < r < R = 2, and we use the space mesh size Ar = 0.035 and a CFL number equal 
to 0.9. The time-asymptotic solution correspond to the data pR = 528.0 and vr = 0.0011. 

The numerical results are plotted in Figure 5.3. Observe that the initial discontinuous pertur- 
bation evolves, gets smoothed out, and spreads in both directions. In particular, in Figure 5.3 (a) 
and (b), we compare the numerical solutions for very large times with the initial steady solution, 
as well as with the steady solution determined from the values V{R) and p{R) of the numerical 
solution. Hence, we have demonstrated numerically that solutions converge to steady solutions 
and, in the present test, the initial perturbation has the effect of modifying the steady solution of 
reference. 

We observed numerically that the solution reaches another steady state, which we plot on the 
ssame figure, for the sake of comparison. The late-asymptotic solution is found to be p{R) - 528.0 
and v{R) - 0.0011. Furthermore, in Figure 5.4, we have computed the relative numerical error in 
a log-log scale in terms of the ratio ^. The convergence rate for the first scheme was foimd to be 
ui = -0.83, while fl2 = -1.43 for the second scheme. 

6 Concluding remarks 

We have presented a well-balanced scheme for relativistic hydrod5mamics posed on a fixed 
backgrotmd spacetime and, especially, Schwarzschild spacetime. For simplicity, we assumed 
that the equation of state is given by a linear relationship between the mass-energy density and 
the pressure. The generalization to other pressure laws is possible; it would lead to "highly 
nonlinear" algebraic expressions, but would not bring any new difficulty for the piurpose of this 
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Figure 5.1.2. Asymptotic solution with the well-balanced scheme. 
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(c) (d) 
Figure 5.2.1. Propagation of a shock/rarefaction pattern (two schemes). 



paper. 

To encompass other backgrounds (for instance, Schwarzschild-de Sitter spacetime, with cos- 
mological constant included), the analysis in Section 3 should be revisited. An analogue charac- 
terization of steady solutions could be derived without additional conceptual difficulties, so that 
our method appears to be relevant for a class of background black hole spacetimes. 
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Although the proposed finite volume method is currently restricted to problems involving 
one spatial variable, it is of genuine interest for numerical relativity and provides a new tool in 
order to investigate the evolution of self-interacting matter under symmetry conditions. Indeed, 
a suitable extension of our method IITOl allows one to deal with the coupled Einstein-Euler system, 
in which the metric itself is an unknown of the problem. This strategy has now been applied 
to the Einstein equations for spherically symmetric spacetimes and should also be useful to 
investigate T^-symmetric matter spacetimes (admitting, by definition, two commuting Killing 
fields). Specifically, we refer to |10| for a study of the spherically symmetric, self-gravitational 
collapse of compressible fluids (investigated earlier by Novak and Ibanez [17J and Papadopoulos 
and Font ||T9| in a different gauge) and for the application of the proposed finite volume technique 
to the formation of trapped surfaces. 

The treatment of the full Einstein system without symmetry assumptions is currently out 
of the scope of the existing techniques. An intermediate goal will be to encompass general 
fluid equations (without symmetry and with general pressure laws) and arbitrarily curved four- 
dimensional background geometries (satisfying Einstein equations). 
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